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Abstract 

We consider a model of a dilute Bose-Einstein condensed gas at finite tem- 
peratures, where the condensate coexists in a trap with a cloud of thermal 
excitations. Within the ZGN formalism, the dynamics of the condensate is 
described by a generalized Gross-Pitaevskii equation, while the thermal cloud 
is represented by a semiclassical kinetic equation. Our numerical approach 
simulates the kinetic equation using a cloud of representative test particles, 
while collisions are treated by means of a Monte Carlo sampling technique. 
A full 3D split-operator Fast Fourier Transform method is used to evolve the 
condensate wavefunction. We give details regarding the numerical methods 
used and discuss simulations carried out to test the accuracy of the numerics. 
We use this scheme to simulate the monopole mode in a spherical trap. The 
dynamical coupling between the condensate and thermal cloud is responsible 
for frequency shifts and damping of the condensate collective mode. We com- 
pare our results to previous theoretical approaches, not only to confirm the 
reliability of our numerical scheme, but also to check the validity of approxi- 
mations which have been used in the past. 

PACS numbers: 03.75.Fi, 05.30.Jp, 02.70.Ns, 67.40.Db 
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I. INTRODUCTION 



Bose-Einstein condensation (BEC), whereby bosons form a condensate by macroscopi- 
cally occupying the lowest energy state of the system, is a striking and important consequence 
of quantum statistics at low temperatures. The resultant long-range order manifests itself 
in phenomena such as macroscopic coherence and superfluidity. In general the condensate 
is depleted by correlation effects and through thermal population of excited states at finite 
temperatures. The former, termed quantum depletion, is particularly important for dense 
fluids such as liquid 4 He, where only around 10% of the atoms are condensed in the low 
temperature limit. In contrast, the quantum depletion in trapped, dilute gaseous BECs 
is typically less than 1% [d,^]. The noncondensed fraction is thus mainly composed of 
thermal excitations, and almost pure condensates can be prepared by evaporative cooling 
to very low temperatures. Atomic vapors therefore allow unique opportunities to study the 
properties of Bose condensates under a wide range of conditions, from the pure condensate 
phase to the noncondensed thermal cloud above the BEC transition. 

The condensate in a dilute Bose gas is well-described by means of a macroscopic wave- 
function, which in the limit of low temperatures evolves according to the Gross-Pitaevskii 
(GP) equation. Well-known techniques allow both numerical and analytical solutions of this 
equation, and comparisons with experiment at low temperatures show excellent agreement 
for both static and dynamical properties [f|]. However, generalizations of the theory to finite 
temperatures, where thermal excitations coexist with the condensate, have proved far more 
difficult. To accurately describe the dynamical behavior in this situation requires a theory 
which treats both components in a fully consistent manner. Such theories have recently 
been formulated, but the challenge of obtaining explicit solutions has remained. What has 
been lacking in particular is a computationally feasible method for treating the dynamics of 
the thermal cloud. It is these computational aspects that concern us most in this paper. 

The earliest studies of dynamics at finite temperatures were based on the Hartree-Fock- 
Bogoliubov (HFB) approximation [^,0. Within this theory, excitations of the condensate 
are obtained by solving the HFB equations which are derived by linearizing the GP equation 
about the equilibrium solution, or equivalently, from the grand canonical Hamiltonian of the 
system JF| . The frequencies of the excitations are identified with the collective modes of the 
condensate. This theory, however, is incomplete. Although the excitations are thermally 
populated, the condensate in fact oscillates in the presence of a static thermal cloud. This 
ignores the dynamical response of the thermal cloud to condensate fluctuations which is 
responsible for Landau damping and associated frequency shifts. By the same token, the 
theory cannot be used to account for the response of the system to external perturbations 
as typically used in experiments to excite the trapped gas P~PT]. This problem becomes 
critical at high temperatures, where collective motion of the thermal cloud can exert a 
major influence on the condensate evolution, as reflected in experimental results for the 
mode frequency and damping rate. 

Recent important work by Morgan et al. [|T2|,r3| and Giorgini [13,13] has extended the 
HFB theory to include collisionless noncondensate dynamics within second-order perturba- 
tion schemes, and derived expressions for damping rates and frequency shifts of low-energy 
modes. A variant of these approaches is the dielectric response formulation of Reidl et al. 
T6fl . One limitation of these theories is the absence of collisions which require a kinetic 
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theory for their description. Quantum kinetic equations for BECs have been developed by 



Gardiner and collaborators 1171, Stoof 181 and Walser et al. 



However, calculations 

based on these theories are very difficult to carry out and as a result, they have not yet 
been used to study collective excitations. A somewhat simpler scheme is the one developed 
by Zaremba, Griffin and Nikuni (ZGN) pU,2"I|, which treats the excitations semiclassically 



within the Hartree-Fock (HF) and Popov approximations. One can then identify the excita- 
tions with a thermal cloud of particles, with dynamics governed by a Boltzmann equation for 
the phase-space density. In analogy with its classical counterpart, binary collisions between 
particles are described by means of a collision integral; however, an additional collision in- 
tegral arises to account for collisions with the condensate. The latter leads to an important 
modification of the GP equation which must now include a non-Hermitian source term to ac- 
count for the transfer of atoms into and out of the condensate. This process, taken together 
with mean field coupling between the two components, leads to damping and frequency shift 
of the condensate collective modes at finite temperature. 

The coupled GP and Boltzmann equations are far from trivial to solve, and several 
approximations have been invoked in the literature in order to explore their properties. When 
the characteristic collisional time scale, r, satisfies uqt <C 1, where uq is the trap frequency, 
then collisions dominate and the system is said to be in the hydrodynamical regime. One 
can then take moments of the kinetic equation to derive a set of coupled hydrodynamic 
equations for the noncondensate which, together with similar equations for the condensate, 
can be solved under certain conditions [2l]J22|. In the opposite collisionless regime, ujqt 3> 1, 



Stoof and co-workers |23]j24| have used a joint variational and moment scheme to model the 
condensate and noncondensate, respectively, while Nikuni [25] recently applied a moment 

,[2j],[l(J. Although these moment methods provide 



method to study the scissors mode 



some insight into the coupled dynamics of the two components, they constitute a truncated 
description which precludes coupling to internal degrees of freedom of the gas. Thus, Landau 
damping is neglected. In order to avoid this limitation, and to facilitate direct comparisons 
with experiment, one must resort to the full kinetic theory. It is therefore desirable to 
directly simulate the ZGN equations without making approximations beyond those used to 
derive the equations themselves. In this paper, we describe a technique to calculate the 
dynamics of the thermal cloud using iV-body simulations. Within this approach, a swarm 
of test particles is used to represent the evolution of the semiclassical phase-space density, 
while collisions are handled using a Monte Carlo sampling technique. The dynamics of the 
condensate on the other hand is determined by numerically propagating the GP equation 
using a split-operator Fast Fourier Transform (FFT) method. Application of the method 
to the quadrupole and scissors [ PU| , PU[ modes has been discussed elsewhere, and in 



both cases, good agreement with experiment ||8|JT0|] was found. Although an outline of the 
numerical methods used was given in this earlier work, we give much more detail in the 
present paper. 

This paper is organized as follows. In Sec. |IJ we briefly review the ZGN formalism, before 
discussing our numerical methods in Sec. P|. In Sec. [TV] the Monte Carlo sampling is tested 
by comparison of equilibrium collision rates against semi-analytic calculations. Landau and 
collisional damping rates for the monopole modes in spherical traps are also compared to 
previous theoretical treatments. We sum up and outline possible future research directions 
in the Conclusion. 
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II. THE ZGN FORMALISM 



We begin by reviewing the ZGN formalism, which was derived and discussed in detail in 
Ref. For a Bose-condensed gas one can decompose the second-quantized field operator 
ip(r,t) in the following manner 

^(r,t) = $(r,t) + ^(r,t), (1) 

where the ensemble average $(r, t) = (^(r, t)) takes on a non-zero value due to Bose broken 
symmetry, and is identified with the condensate wavefunction. The remaining field operator 
i/j(r, t) has a zero expectation value and corresponds to the noncondensed component of the 
cloud. The second-quantized Hamiltonian for the system is given by 

H= /dr^(r) -^- + U ett (r) ^(r) 

+ l -J drdr'$\r)^\r')U int (r,r')^(r')i>(r), (2) 

where in most cases the trap is well approximated by a harmonic potential U cxt (r) = 
m{oj 2 x 2 + oj 2 y 2 + uj 2 z 2 )/2. We also assume a contact interaction, U int (r,r') = g<5(r — r'), with 
g = 4irh 2 a/m, where a is the s-wave scattering length and m is the atomic mass. Using 
ihd t ij) = [ip, H] with ([!]) and (Q), one can derive coupled equations of motion for the conden- 
sate and thermal cloud. In particular, the condensate order parameter evolves according to 
a generalized form of the GP equation 

d / h 2 V 2 \ 

z£— $(r,t) = \ -^— + U ext (r)+g[n c (r,t) + 2n(r,t)]-iR(r,t)J^r,t), (3) 

where n c (r,t) = |$(r,t)| 2 and n(r,t) = (^(r, t)i()(r, t)) are the condensate and nonconden- 
sate densities respectively. In arriving at this equation we make the Popov approximation 
whereby the so-called "anomalous" density, m(r, t) = (ip(r,t)ip(r,t)), is neglected. This 
sidesteps problems associated with including this term, such as ultraviolet divergences and 
an unphysical gap in the energy spectrum at low momenta [0. To go beyond this approxima- 



tion in a consistent manner requires a careful treatment of interparticle collisions [jrj| , and is 
beyond the scope of the present work. The source term R(r, t) is an important modification 
of the usual GP equation as it allows the normalization of the wavefunction $ to change 
with time. Physically this is due to collisions between condensate and noncondensate atoms 
which have the effect of transferring atoms into or out of the condensate. The source term 
will be defined in terms of a collision integral later. 

It is convenient to describe the dynamics of the noncondensate in terms of the Wigner 



operator pT|j31| , which leads to the definition of a phase-space distribution, /(p, r, t), for 
the thermal excitations. The equation of motion for the noncondensate can then be written 
as a kinetic equation 

-/(p,r,t) + £ ■ V/(p,r,t) - W(r,t) ■ V p /(p,r,t) = C 12 [f] + C 22 [/]. (4) 
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In deriving this equation a number of approximations have been made, some of which have 
already been mentioned. Importantly, the excitations are assumed to be semiclassical within 
the HF approximation; an excitation with momentum p possesses an energy e = p 2 /2m + 
U(r,t), where the effective potential U(r,t) = U ext (r) + 2g[n c (r,t) + h(r, t)} is composed of 
the trap potential as well as mean fields from the condensate and the thermal cloud. The 
noncondensate density appearing in this expression is given in terms of the distribution by 

S(r - () = /(2^ /(P ' r - t) - (5) 

The terms on the right-hand side of (^) are collision integrals that represent binary 
collisions between atoms. The C 22 term is familiar from the kinetic theory of a normal Bose 
gas, and corresponds to the scattering of two atoms from initial to final thermal states. It 
is given by 

C 22 [/] = , o , / dp 2 dp 3 dp 4 5(p + p 2 - P3 - P4) 

nh 6 m z J 

x5(e + e 2 - e 3 - e 4 )[(l + + / 2 )/ 3 / 4 - // 2 (1 + / 3 )(1 + / 4 )]. (6) 

where / = /(p, r, t) and fa = f(pi, r, t). The total bosonic cross-section is given by a = 8na 2 . 
The delta functions enforce momentum and energy conservation in the collision, while the 
factors (l + /j) account for Bose-enhancement of the scattering. The analogous Ci 2 collision 
integral corresponds to collisions that involve a condensate atom in either the initial or final 
states. It is given by 

Cn[f] = ——5 / dp 2 dp 3 dp 4 5(mv c + p 2 - p 3 - p 4 )5(e c + e 2 - e 3 - e 4 ) 
Tcm z J 

x [S(p - pa) - 5(p - pa) - S(p - p 4 )][(l + f 2 )f 3 f, - / 2 (1 + / a )(l + fa)}, (?) 
where the local condensate velocity and energy are respectively given by 

2tm\9r 



and 



1 2 

-mv c + fi c . 



Here, \x c is the condensate chemical potential defined as 

h 2 v 2 jw c TT 

He = + U cxt + gn c + 2gn. 

2m y/n c 

If the condensed and noncondensed components are in local equilibrium, the Ci 2 integral 
vanishes. Conversely, when the system is perturbed from equilibrium the C\2 term acts to 
transfer atoms between the condensate and thermal cloud. These collisions then define the 
source term in (|3p according to 

The relative numbers of condensate and thermal particles will then adjust as a function of 
time until local equilibrium is re-established. 
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III. NUMERICAL METHODS 



In this section we describe the numerical methods used to solve the ZGN equations (|3[)- 
(§) in the context of a dynamical simulation. We first discuss the numerical methods used 
to solve the GP and collisionless Boltzmann equations. Although these are based on well- 
established techniques (see e.g. we feel that our partly pedagogical discussion will be 
useful for those trying to reproduce our simulations, while highlighting the correspondence 
between the quantum and classical dynamics of the system. We then move on to discuss 
treatment of the C22 and Cyi collision integrals by Monte Carlo sampling. Finally, an 
overview of the simulations is provided, including a discussion of how one calculates the 
equilibrium initial state of the system, as well as estimating the phase-space density in real 
time for use in evaluating the collision integrals and (|7|). 



A. The Gross-Pitaevskii equation 

For the benefit of the following discussion we rewrite the GP equation (§) in the form 

ih^<f>(t) = H(t)Q(t) . (9) 

The time dependence of the Hamiltonian, H(t) = T + V(t), arises from the potential V(t) 
which also includes the non-Hermitian source term R(r,t). In most of our simulations the 
time dependence is dominated by the nonlinear condensate potential and it is this term 
which is the main source of numerical instabilities when the number of condensate atoms is 
large. It is therefore important to develop a numerical algorithm which is accurate even in 
this limit, and at the same time, numerically efficient. 
A formal solution of the above equation is given by 

§(t + At) = U(t + At,t)$(t) (10) 
where the evolution operator U has the expansion 

1 rt+At 1 rt+At rt' 

U(t + At, t) = 1 + — / dt'H{t') + — ^ / / dt'dt"H{t')H{t") + ■■■ . (11) 
Expanding the Hamiltonian as a Taylor series, 

"<o = *w + f ('-«) + ^c- «)' + •■■ 

= a + / j(f_ t ) + i 7 (f_ t )2..., (12) 



we obtain 



U(t + At, t) = 1 + ^At + £-(Atf - ^-(Atf + 0{At% (13) 



The lowest order exponential approximant to this expansion is 
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U(t + At,t) ~ e" 



-iH(t)At/h 1 ^£ f\j.\2 

2h dt [ ' 



0(At 3 



(14) 



The error of second order is shown explicitly. The first term on the right hand side is of course 
exact for a time independent Hamiltonian but significant errors arise when the Hamiltonian 
is time dependent. These errors can be minimized by reducing the time step At, but at the 
expense of increasing the computation time required to complete a simulation. Since this 
imposes practical limits on the physical problems that can be addressed, a more accurate 
approximant is desirable. 

A higher order exponential approximant is provided by 



U(t + At,t) ~ e 



-i(a+±/3At)At/h 



(15) 



A comparison with ( [TBI ) indeed confirms that the error is 0(At 3 ). To this order of accuracy, 
we can make use of ([T2|) to estimate (3 by reverse differencing, 



H(t) - H(t - At) 
At 



and thus obtain 



where 



U(t + At, t) ~ e -^(t)At/h + ( At 3 



(16) 



(17) 



H{t) =T + V(t) 



with 



V(t) 



3V(t) - V(t - At) 



(18) 



(19) 



This is recognized as an approximation to the potential at time t + At/2, the midpoint of the 
current time step, as obtained by a linear extrapolation from the potential at times t — At 
and t. 

The implementation of (17-19) is very simple and costs only a small additional amount of 
memory to store the potential from the previous time step. The actual numerical represen- 
tation of the evolution operator can be achieved by various methods. One popular approach 
is the Crank-Nicholson method |R3 , where finite-differencing Cayley's form for the operator 
leads to a set of linear equations for the wavefunction at discrete grid-points in r. The 
problem then reduces to decomposition of a tridiagonal matrix at each time step and along 
each spatial dimension. In contrast, we favor a split-operator method, where a factorization 
of the exponential is effected by means of the Baker-Campbell-Hausdorff (BCH) formula. 
One finds that 



-iHAt/K 



e -iVAt/2h Q -iTAt/h e -iVAt/2h _|_ 



AtY 
12 \~ih) 



T, V 



T + 



V 



+ C(At 4 ). (20) 



The error generated by this approximation is of the same order as found in (0). Applying 
the first term on the right hand side then evolves the wavefunction to second order accuracy 
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in At. In principle, higher order schemes can be constructed by splitting into more elab- 
orate combinations of the V and T operators. However, to justify the effort, an improved 
approximation for H(t) is required. We have found that second order accuracy is sufficient 
for most applications, although difficulties do arise if the time scale of the simulations is 
exceedingly long. 

The split-operator scheme fl20|) is straightforward to implement with a discrete grid 
in position space. The two potential steps are applied by multiplying the wavefunction 
at each grid point by e - lVAt / 2h ^ while the kinetic term e - lTAt / h [ s conveniently treated in 
momentum space. The limiting step in the calculation is therefore the application of forward 
and inverse Fast Fourier Transforms (FFTs) at each time step, but efficient FFT routines 



for arbitrary numbers of dimensions are readily available |35fl . The dynamical evolution of 
the wavefunction can thus be followed over a series of time-steps. Alternatively, stationary 
solutions of the time independent GP equation can be easily found by evolving the time- 
dependent equation in imaginary time t — > —it. 

A typical application provides some indication of the relative merits of the higher order 



approximant in ( [T7| ) as opposed to the lower order scheme in (|H|). With the latter, one finds 
a monotonic increase in the energy expectation value with time. In simulations of a collective 
mode this effect would be apparent as a slow increase in the mode amplitude, which is clearly 
undesirable when quantifying damping at finite temperatures. More importantly, since the 
rate of increase scales with the mode energy, higher frequency excitations tend to build 
in amplitude more rapidly. These excitations are initially generated at a low level by the 
numerics; however over sufficiently long simulation times they eventually lead to instabilities 
in the wavefunction. These problems are essentially eliminated with the higher order scheme. 
The stability of the simulations is dramatically improved and the energy tends to oscillate 
with small amplitude about a constant value, rather than increasing monotonically. The 
improved stability allows much larger time-steps to be taken without compromising accuracy, 
leading to a considerable saving in computational effort. 



B. Collisionless particle evolution 

In this section we discuss solution of the collisionless Boltzmann equation (C12 = C22 — 0) 
using iV-body simulations. The effect of collisions is dealt with later. Collisionless Boltz- 
mann (or Vlasov) equations which include mean-field interactions arise in many disparate 
fields, such as plasma physics, condensed matter physics and astrophysics. Since the equa- 
tion involves phase space variables in six dimensions, it is generally very difficult to solve 
using standard methods for treating partial differential equations. An alternative approach 
used extensively in the literature is to represent the phase-space density /(p, r, t) by a cloud 
of discrete test particles . The momentum and position of each particle in an external po- 



tential U(r, t) is then evolved according to Newton's equations. The phase space distribution 
for this situation is given by 

1 1=1 

where the weighting factor is fixed by the requirement that the phase-space distribution is 
normalized to the number of physical atoms, N, with Nh 3 = J drdp/. By using a sufficiently 
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large number of test particles, Nt, a reasonable approximation to the continuous phase space 
distribution is obtained. Note that the number of test and physical particles is not necessarily 
equal. In fact, for a relatively small number of physical atoms (N ~ 10 4 ) it is essential to 
simulate more test particles (Nt > 10 5 ) in order to minimize the effects of a discrete particle 
description. Conversely for large samples one can simulate fewer "superparticles" so that 
the calculations are not too intensive. 

The phase space variables are updated by advancing the position and momentum of each 
particle at discrete time steps At. This is not as trivial as one might naively expect. Conven- 
tional integration schemes for ordinary differential equations, such as classical Runge-Kutta 
methods, can lead to non- conservation of energy over long-time simulations when applied 
to Hamiltonian systems. This results in spurious damping or excitation of the system. In 
contrast, symplectic integrators [0,0 are use d extensively in molecular dynamics (MD) 



simulations since they possess several desirable properties, such as conservation of phase- 
space volume and of energy over a long period (as is required in autonomous Hamiltonian 
systems). We use a second-order symplectic integrator in our calculations, which is the 
classical analogue of the split-operator method discussed earlier. To show this, it is conve- 



nient to work within the Lie formalism |33| . Consider the classical Hamiltonian for a single 
particle, Hi = pf/2m + V(r»). The evolution of its phase-space coordinates Zj = (pi,r$) is 
then determined by the equation 

^ = K H t } = -iCz h (22) 

where {F, G} = J2j d rj Fd Pj G — d Pj Fd rj G is the Poisson bracket and C is the Liouville 
operator ||38|| . One can then write 

z(t + At) =e~ iCAt z(t). (23) 

Splitting the Hamiltonian into potential and kinetic terms, Hi = T(pi) + V"(r$), the BCH 
formula can be used again to show that [37] 



One now sees the analogy with the quantum operator (0), where both conserve energy to 
order (At) 2 . The effect of the classical operator (pi) in the simulations is to update the 
particle positions and velocities in three steps 

f< = Ti(t) + iAtVi(t), 

Vi(t + At) = Vi(t) - m' 1 AfVy(f ! ), 
Ti(t + At) = ii + |At Vi (t + At). (25) 

By analogy with (|I"9l), V should be the midpoint value of the potential, V(t), when it is 
time-dependent. In our simulations, V is the effective potential U(r,t) = t/ ext (r) + 2gn(r,t) 
felt by the thermal atoms, where n = n c + h is the total density. 
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C. Thermal Cloud Potential 



The effective potential U is determined self-consistently as the system evolves in time, and 
includes the condensate mean field 2gn c (r, t) and the mean field generated by the thermal 
cloud 2gh(r, t). The latter is in general much weaker than the condensate mean field due to 
the larger spatial extent (and therefore lower density) of the thermal cloud. Nevertheless, it 
is important to include this term in order to ensure the conservation of the total energy of 
the system. In addition, from the perspective of the condensate, the noncondensate mean 
field is necessary in order to account for the temperature-dependent damping and frequency 
shifts of condensate collective modes. 

Although the calculation of the condensate mean field is straightforward, the use of 
discrete particles with a contact interatomic potential creates a problem in determining the 
noncondensate mean field. Taken literally, the mean field consists of a series of delta peaks 



This expression clearly cannot be used as it is to generate the forces acting on the test par- 
ticles that are required in the MD simulation. Rather, the density h T (r,t) must be replaced 
by a smooth and differentiable thermal cloud density and some smoothening operation is 
therefore needed. A possible first step might be to divide space into cells and to determine 
the mean density within each cell by binning the test particles appropriately. However, this 
binning procedure generates spatial discontinuities on the scale of the 3D grid being used 
that would still have to be smoothed out in some way. In addition, temporal discontinu- 
ities arise as particles migrate from one cell to another. These temporal fluctuations are of 
course spurious since they depend on the number of test particles and decrease in relative 
amplitude as this number is increased. It is apparent that the binned density has some un- 
desirable properties associated with the statistical fluctuations in the number and positions 
of particles in each cell. 

As an alternative to this binning procedure, we generate a smooth thermal cloud den- 
sity by performing a convolution with a sampling (or smoothening) function S(r) which is 
normalized to unity. In particular, we define 



where we choose S(r) ~ e - r2 h 2 ; i. e . an isotropic Gaussian sampling function of width 77. 
Since V>S| r= o, no force is exerted by a particle on itself and the sum can extend over all 
particles in the ensemble. Ideally, the width of S*(r) should be small compared to the 
curvature of the noncondensate density. If, at the same time, the number of particles 
contributing to the sum at a given position r is large, it is clear that the sampled potential 
will be relatively smooth. Note that the smoothening operation is equivalent to assuming a 
finite-ranged interatomic potential. 

The sampled potential (or its gradient) is needed at the position of each test particle 
and at the mesh points on which the condensate wavefunction is defined. However a direct 




(26) 
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summation for all points would be prohibitive. We therefore proceed by making use of a 
FFT. First, each particle in the ensemble is assigned to points on the 3D Cartesian grid 
using a cloud-in-cell method |!6|. This is most readily explained in ID: consider a particle 
at position x, between two grid points at x/. and Xk+i- The particle is assigned to both 
points with weightings (1 — a) and a respectively, where a = (x — Xk) / (xk+i — x k)- This can 
be viewed as a more sophisticated binning procedure in that it takes into account the actual 
positions of particles within the cells. The generalization to 3D is straightforward, where in 
this case the particle is assigned to the eight points which define the unit cell containing the 
particle. We then convolve the cloud-in-cell density with the sampling function by Fourier 
transforming it and then multiplying it by the analytic FT of the sampling function. An 
inverse FFT then generates the sampled potential on the 3D grid. This potential is used 
directly in the GP evolution, while the forces on the test particles are obtained by taking a 
numerical derivative and interpolating to the positions of the particles. 

This overall scheme is illustrated in Fig. |I[ The solid line shows the equilibrium thermal 
cloud density along a line through the center of an isotropic trap with trapping frequency 
ujq = 2n x 187 Hz, a system we study in more detail later. The trap contains a total of 
A/tot = 5 x 10 4 87 Rb atoms, and at a temperature of T = 250 nK there are JV ~ 4.0 x 10 4 
thermal atoms. The rapidly fluctuating dashed line is the density along this line produced 
by the cloud-in-cell method using a thermal distribution of iVV — 4.0 x 10 5 test particles, 
that is, ten times the actual number of thermal atoms. The effect of statistical fluctuations 
is clearly evident. Finally, the smooth dashed line is the result of the convolution using a 
width parameter of rj ~ 0.76 a^, where a\ lo = (h/mtoo) 1 ^ 2 ~ 7.9 x 10 _7 m is the harmonic 
oscillator length for the trap being considered. (For comparison, the mesh size is Ax ~ 
0.27ah o .) It should be noted that the dramatic smoothening of the density achieved is 
partly a consequence of performing a full 3D convolution; a ID convolution of the cloud-in- 
cell density with the same width parameter would not reduce the amplitude of the spatial 
fluctuations to the same degree. Finally, we compare the convolved density to the actual 
equilibrium density. Apart from differences due to the statistical sampling of test particles, 
one can see that the peaks in the thermal cloud density at the edges of the condensate are 
slightly broader in the convolved density, as would be expected. However, the differences 
are minor and do not affect the dynamics of the system significantly. For consistency, the 
n c term appearing in U(r, t) is also convolved. 



D. Collisions 

The methods outlined so far allow one to follow the condensate wavefunction and tra- 
jectories of the atoms subject to a time-dependent potential, so long as the system is in 
the collisionless regime. However in general the collisional terms in the Boltzmann equation 
will be nonzero C22 7^ 0, G\i 7^ 0. In other words, during each time step there is a certain 
probability that a given test particle will collide with another thermal atom or with the 
condensate. If the typical collision timescale r is such that r ^> At, one can treat the free 
particle evolution and collisions separately. Each particle's trajectory is first followed using 
the methods discussed in the previous section, and the possibility of collisions occurring is 
then considered at the end of the time step. Probabilities for either C22 or C\2 collisions 
are calculated in a way which is consistent with a Monte Carlo sampling of the collision 
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integrals, as discussed below. 



1. C22 collisions 

We first give details for the C22 integral (|B|), which physically corresponds to scattering of 
two thermal particles into two final thermal states. Hence the process conserves the number 
of thermal atoms, J dp/(27r/i) 3 C22 = 0. We are interested in the mean collision rate at a 
point r (as defined in Appendix A), which is given by 

II? = n f^ m 2 J j d P 2 ^ 2 J d P3 j d P4^(Pl + P2 - P3 - P4) 

x5( ei + e 2 - e 3 - e 4 )(l + / 3 )(1 + h)- (27) 

For our purposes it is convenient to express the integral in terms of new momentum variables 
(p , p'o) and (p', p"): Pi, 2 = (po ± p')/V2 and p 3j4 = (p' ± p")/V2- Po and p' are 
proportional to the center-of-mass and relative momenta, respectively, of the incoming 1 
and 2 particles. By implicity assuming energy and momentum conservation (po = p(j, 



p' = p") one can rewrite (27) in the simplified form 



rif = / [ T^U I ^-Ivi - v 2 |(l + /a)(l + h) , (28) 



(2ttH) 3J J (2irh) 3J J 4vr 

where p 3j4 = (po ±p'u(fi))/v / 2, with u(f2) a unit vector in a direction specified by the solid 
angle Q. Calculation of the rate therefore involves integrals over all possible initial states 
and all scattering angles Q. In the equilibrium situation, this rate defines a local mean 
collision time r° 2 according to 



n 



^22 — ' (29) 
r 22 



where no(r) is the equilibrium thermal cloud density. As shown in |3(J, l/r^ below T c 
is a strong function of position for a trapped Bose gas and is peaked at the edge of the 
condensate. In the classical (i.e. Maxwell-Boltzmann) limit, 1/t°2 reduces to v^crw^fio, 
with v th = (SkT/irm) 1 / 2 . 

To relate this to collision probabilities for individual atoms in our simulations requires 
sampling of the integral using a rejection method as discussed in detail in Appendix A 
|34] , |3"9"|| . At each time step atoms are first binned into cells of volume A 3 r according to their 
position. The atoms within each cell are then paired at random, and a probability for a pair 
(ij) to collide in the time step At is assigned according to 

P» = fi ff \y i -y j \ J ^(l + / 8 )(l + / 4 )Af. (30) 

The integral over Q can be evaluated by averaging over a sample of randomly selected final 
states which are obtained by choosing uniformly-distributed random values for the scattering 
variables cos# and <fi. However, in simulating the collision process, the velocities of the 
incoming particles must actually change to a specific, but random, pair of final velocities. 
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These velocities lie on a sphere centered at (vi + v 2 )/2 with a radius |vi — v 2 |/2 and can be 
chosen by randomly selecting the scattering angle Qr. The appropriate collision probability 
for this event is then 

P 22 = Hv* - v,|(l + f^)(l + f?«)At. (31) 

This probability depends upon the phase space densities of the final states, f^ R , f^ R } reflect- 
ing Bose statistics. If this single scattering probability is averaged over a random distribution 
of scattering angles Qr we recover the average probability defined in (^) . 

The simulation of C 22 collisions thus proceeds as follows. A pair of test particles (ij) 
in a given cell is chosen at random. Whether a collision of this pair occurs is then tested 
by comparing P? 2 to a random number X 22 uniformly distributed between and 1. If 
X 22 < P 22 the collision is accepted, and the velocities of the test particles are updated 
accordingly. If X 22 > P? 2 , no collision occurs and the velocities of the colliding pair are 
unchanged. In either case, another pair is randomly selected and the procedure is repeated 
for all pairs in each cell of the sample. 



2. C\2 collisions 

The C12 collisions are treated in a similar manner to C 22 . The key difference here is that 
one of the collision partners is a condensate atom in a definite state, and it is necessary 
to distinguish the collisional processes which either transfer an atom into or out of the 
condensate. For example, the "out" collision rate as defined in ([All ) is given by 



irm 2 h 3 



CTl 

T if = I dp 2 dp3dp 4 5(p c + p 2 - Ps - p 4 )5(e c + e 2 - e 3 - e 4 )/ 2 (l + / 3 )(1 + U). 



(32) 



This represents scattering of a thermal atom from the condensate to produce two thermal 
atoms. The reverse process gives the "in" collision rate defined in ( |A15| ): 

r i2 = / dp 2 dp 3 dp 4 5(p c + p 3 - p 2 - p 4 )6(e c + e 3 - e 2 - e 4 )/ 2 (l + / 3 )/ 4 . (33) 

In obtaining (^) we have interchanged the 2 and 3 labels in order to define an integral 



having the same / 2 weighting factor as in (|32|). These two integrals give the true "in" and 



"out" collision rates. However in the simulations, it is useful to drop the cubic terms / 2 / 3 / 4 
which formally cancel exactly between the "in" and "out" rates. Since these two rates are 
evaluated differently as explained below, this cancellation will not be numerically precise, 
and it is therefore preferable to eliminate the cubic terms from the calculation of collision 
probabilities. In the following, we denote the rates with the cubic terms removed by ° ut . 
Dropping these terms of course does not change the net rate of transfer from the condensate 
to the thermal cloud that actually takes place. 

The "out" term can be reduced by transforming the momentum variables as before, with 
the result 
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r " = i0f bn ^S^ (i+h+h) ' (34) 



where v° ut = y |v c — V2I 2 — 4gn c /m is the relative velocity of the initial states, corrected to 
account for energy conservation (locally, the mean field energy of a thermal atom is higher 
than that of a condensate atom by an amount gn c ). Now, if we consider each atom in the 
distribution / 2 in turn, the probability for collision with the condensate is given by 



= n c av™\l + ff R + ft*) At . (35) 

In this case, the final thermal atom velocities V3, V4 lie on a sphere of radius t>° ut /2 centered 
on (v c + v 2 )/2, with a random scattering angle Qr. 

"In" collisions involve scattering of two thermal atoms to produce a condensate and a 



thermal atom. In the context of (|33|), the incoming atoms are labelled 2 and 4, and the 
outgoing thermal atom is labelled 3. Energy-momentum conservation in (|33"D dictates the 
condition (p c — p 2 ) • (p c — p 4 ) = mgn c . Thus, unlike the case of C 22 collisions, one cannot 
arbitrarily select a pair of 2 and 4 atoms from the sample since this condition will in general 
be violated and the collision cannot occur. To proceed, we perform the integrations involving 
the delta functions in (|33"D to obtain 

15 - f M~J^ I d*A. (36) 



{2nhf J -nv\, 

where v* n = v 2 — v c is the velocity of thermal atom 2 relative to the local condensate velocity. 
The second integral is a two-dimensional integral over a velocity vector v which is in a plane 
normal to v™. The velocity of the other incoming thermal atom, particle 4, is given by 

, - , 9n c .in 
v 4 = v c + v + -v r , 

while the velocity of the outgoing thermal atom is 

g^c .in 



V 3 = V 2 + V + 

In the simulation one considers each thermal atom in the distribution / 2 in turn, then 
randomly selects two numbers that define the vector v = vr within a plane of area A v . The 
collision probability is then given by 

Pr = ^fZ«At. (37) 

Note that the area A v appears in this expression, which at first sight is disconcerting since it 
is an arbitrary number entering as a simulation parameter. However, we find that the total 
rate is largely independent of this area so long as the plane completely samples the occupied 
regions of phase space. We show results confirming this statement in the following section. 

This analysis yields probabilities for a particular atom to undergo "out" or "in" collisions. 
To decide whether either event takes place, another random number < X 12 < 1 is chosen. 
If X 12 < P° ut then an "out" collision is accepted; the incoming thermal atom is removed 
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from the ensemble of test particles and two new thermal atoms are created. On the other 
hand, if P° ut < X 12 < P° ut + P- n , then an "in" collision takes place and atom 2 is removed 
from the thermal sample. In addition, a second test particle, atom 4, is removed and a 
new thermal atom, atom 3, is created. In practice, it is exceedingly unlikely that a test 
particle will exist that will precisely match the required phase-space coordinates of particle 
4. We therefore search for a test particle in neighboring phase-space cells and remove this 
particle if one is found. This can be justified by remembering that we are only interested 
in describing the evolution in phase space in a statistical way — it is misleading to think of 
a direct correspondence between the test particles and physical atoms. If no test particle 
exists in the vicinity of v 4 , the local phase-space density / 4 , and hence P™, will be zero and 
the "in" collision is precluded from occuring in any case. 

The above procedure leads to a change in the number of atoms in the thermal cloud. 
In order to conserve the total particle number the GP equation (Q) is propagated with 
the i?-term which changes the normalization of the wavefunction and hence the condensate 
number. This quantity can be evaluated from the Monte Carlo process decribed above by 
summing probabilities for particles around each grid point using (H), i.e., 



In practice this assignment to grid points is performed with a cloud-in-cell approach similar 
to the one described earlier. Of course the normalization of the condensate wavefunction 
varies continuously as opposed to the variation of the thermal atom number which changes 
by discrete jumps. Nevertheless, one can show that the subsequent change in the condensate 
normalization is consistent with the addition or removal of atoms from the thermal cloud, so 
that the total particle number, N tot , is conserved within statistical fluctuations (~ \A/V tot ). 



So far we have described various aspects of the numerical scheme. The aim of this 
subsection is to tie these disparate elements together with an overview of the simulation 
procedure as a whole. One of the main applications of our approach is to the study of 
small amplitude collective oscillations around the equilibrium state. The first requirement 
of such a calculation is therefore the self-consistent determination of the equilibrium thermal 
cloud distribution and condensate wavefunction. Since the thermal excitations are treated 
semiclassically, the thermal cloud is described by the equilibrium Bose distribution 



where z(r) = exp{/3[/i c — i/(r)]} is the local fugacity and (3 = 1/ksT. It is straightforward 
to show that both the Cyi and Cyi collision integrals vanish in this case. The noncondensate 
density profile can be evaluated from (§) and ([39]) to yield 




(38) 



E. Overview 




(39) 




(40) 
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where A = {2Tth 2 /mk,BT) 1 / 2 is the thermal de Broglie wavelength. The equilibrium conden- 
sate wavefunction is obtained as the stationary solution of (|^), with R = 0, and the cor- 
responding eigenvalue defines the equilibrium chemical potential /i c . Since the condensate 
and thermal cloud are coupled by mean fields, the two components have to be determined 
self-consistently using an iterative procedure. Details of this have been given by several 



authors (see e.g. [z]J or ||40|| ) and will not be repeated here. 

To represent the thermal cloud in the simulations, an ensemble of test particles must be 
defined. In the case of an equilibrium situation, this ensemble should have a phase-space 
distribution which is consistent with the Bose equilibrium distribution in (|39|). This can be 
achieved using the following rejection algorithm |j34| . First, we distribute particles in position 
space according to the density n(r). To do this, we select three random numbers uniformly 
distributed between — r max and r max , defining Cartesian coordinates, r^, of a particle in 
the occupied region of position space. A further uniform deviate is then chosen from R\ G 
[0, n max ], where n max > max{n(r)}, and compared to the density at that point n(rj). If R\ > 
n(i"j), the particle is discarded and another set of position coordinates selected. Otherwise, if 
R\ < n(rj), the particle is accepted and one proceeds to specify its momentum by choosing 
another random number p^ G [0,p max ]. A random number R\ G [0, / max ] (where / max > 
z(vi)/[\ — with z{jCi) the local fugacity) is compared to f(pi,Ti) to decide whether 

the momentum is accepted or rejected. In the case of rejection another pi is chosen, while 
if accepted two random angles are selected 4> G [0, 27r], cos 9 G [—1, 1], which in turn define 
the momentum vector p^. This procedure is repeated until test particles in the ensemble 
are accumulated. Note that we have exploited the spherical symmetry of the equilibrium 
distribution in momentum space. In principle, a similar method can be applied to position 
space if the trap is spherically or cylindrically symmetric. 

A dynamical simulation can be initiated in one of two ways. Either an appropriate 
nonequilibrium initial state is specified, or the system is dynamically excited with the appli- 
cation of an external perturbation. The latter parallels the procedure used experimentally to 
study small amplitude collective excitations, and usually amounts to some parametric ma- 
nipulation of the trapping potential. Although this might be the preferred approach, it is not 
always the most appropriate, especially when the excitation phase requires a prohibitively 
long simulation time. It is then more convenient to impose the perturbation on the initial 
state itself. Here we are guided by the nature and symmetry of the collective mode being 
studied, as well as information gleaned from earlier calculations such as those based on the 
Thomas- Fermi approximation. For example, the nature of the density fluctuation or velocity 
field associated with the mode might be known and it is then advantageous to use this infor- 
mation in defining the initial state. A good example of this is the breathing, or monopole, 
mode in an isotropic trap. In this case the TF mode has a velocity field v = ar. To impose 
this velocity on the condensate one can simply multiply the ground state wavefunction by 
a phase factor exp(imar 2 /2h) . In the case of the thermal cloud, the same velocity field can 
be imposed by adding ar^ to the velocity of the z-th particle in the equilibrium ensemble. 
This procedure will predominantly excite the lowest monopole oscillation. Although higher 
lying modes might also be mixed in to some extent, they have different frequencies and can 
usually be separated from the dominant mode when analyzing the dynamics. 

Returning to the simulation procedure itself, the condensate wavefunction and thermal 
atom phase-space coordinates are updated in each time step At according to the prescription 
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detailed in Sec. [I i [B|. Then, before treating collisions the thermal atoms are assigned to 
cells in position space. These are used for selecting pairs for C22 collisions, as well as being 
further subdivided into momentum space elements in order to estimate the phase space 
density /(p, r) for calculating collision probabilities. Since collisions are treated one cell at 
a time, the phase space density only needs to be calculated and stored for one particular cell. 
The C\2 and C22 collisions are then treated using the Monte Carlo scheme described earlier 
and the momenta and number of thermal atoms (test particles) are updated. Repeating for 
all of the cells yields the quantity R from fl38|) which, when used in the GP propagation 
(Sec. [Ill Afl , continuously evolves the number of atoms in the condensate. For numerical 
accuracy the positional cells should enclose regions of almost constant thermal density and 
fugacity, and are most conveniently treated using a spatial grid which reflects the (elliptical) 
geometry of the cloud. The momentum elements in contrast lie on a Cartesian grid, where 
a cloud-in-cell method allows one to minimize statistical fluctuations while retaining a fine 
grid for precision. 



IV. RESULTS 

A. Equilibrium collision rates 

Our first calculations are not simulations as such, but are instead checks of the Monte 
Carlo sampling technique we use to evaluate the C\ 2 and C 22 collision rates in real time. 
The physical situation we consider corresponds to the one discussed at the end of Sec. |111 C 



namely 5 x 10 4 87 Rb atoms at 250 nK in an isotropic trap. The equilibrium C22 collision rate 
r° 2 can be evaluated numerically directly from the expression in (|2q) using the equilibrium 
distribution function (|39|). The result as a function of the radial coordinate r is shown as 
the solid line in Fig. 2. The equilibrium Cyi collision rates can also be calculated using the 
equilibrium distribution ( |39"D and equilibrium condensate density n c (r) in (^) or (|33|). The 
"in" and "out" rates are in fact equal to each other in equilibrium and will be denoted F^ 
(recall that these rates are calculated ignoring the cubic terms in the full expression). The 
result of the calculation as a function of r is shown as the solid line in Fig. 3. One sees 
that both the C12 and C22 collision rates exhibit a maximum near to the condensate surface, 
where the fugacity z approaches unity and the equilibrium Bose distribution is strongly 
peaked at p = 0. However, in the case of C22 collisions, the tail of the distribution decays 
more slowly since the thermal cloud density extends out to larger radii than the condensate. 

The Monte Carlo calculation of these rates involves a dynamical simulation of a sample 
of test particles moving in the equilibrium effective potential. The collisionless evolution of 
the particles in time provides an ergodic sampling of phase space. At each time step At in 
the evolution, the collision probabilities in (|3T|), ( p5|) and (|37| ) are calculated and summed 
to obtain a realization of the collision rates at a particular instant of time t n . For example, 
for C22 collisions we have 

p22( f \ 

r° (t ) ~ 2 V n — 

where the sum extends over all pairs of test particles in the cell of volume A 3 r. By repeating 
this calculation over M time steps and performing the average 
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1 M 

(■^22) = ~jT7 2^ ^22 (tn) ; 
n=l 

we obtain the Monte Carlo estimate of the collision rate. The same procedure is used for 
the C12 "in" and "out" rates. To obtain histograms of the collision rate as a function of the 
radial position r, we bin the individual collision probabilities according to the positions of 
the colliding pair. The Monte Carlo results presented in Figs. |2] and |^ were obtained with 
only M = 200 time steps of size uioAt = 0.002, which was already sufficient to give good 
statistics. A comparison with the direct numerical calculations shows very good agreement, 
the main error arising from estimating /(p, r, t) in real-time by binning particles into phase- 
space cells. This was confirmed by repeating the simulation but calculating the collision 
probabilities using the actual equilibrium Bose distribution ([39|) rather than the binned 
approximation to it. One can try to improve the binned distribution but there is a trade-off 
between using smaller phase-space cells which would provide a more accurate representation 
of the distribution, and larger cells which contain more particles and thus improve statistics. 
Our choice of cell size tries to optimize these opposing requirements. 

The main observation to be made about Fig. 3 is that the "in" and "out" C12 rates are 
very similar, despite the very different appearance of the probabilities in (|35|) and (|3"?D . Note 
in particular that these results confirm that the "in" rate is independent of the arbitrary 
area A v in (|37|). It is of course important to minimize the difference between these two rates 
since any imbalance implies a net transfer of atoms between the condensate and thermal 
cloud which should not occur in equilibrium. On the other hand, a calculated imbalance 
partly reflects the fact that the equilibrium state we start with is not the "numerical" 
equilibirium state that is consistent with the various numerical approximations being made. 
In fact, we find that when a full simulation is carried out, the system relaxes to a new, slightly 
different equilibrium. In other words, the system automatically adjusts to compensate for the 
numerical approximations. Nevertheless, it is desirable to avoid an imbalance to whatever 
extent possible. Taking the collision rate histograms in Fig. 3 and integrating over r, we 
find a discrepancy between the total "in" and "out" rates of about 1%. This imbalance can 
be minimized by judicious choice of the shape of the phase space elements (see Sec. |Ill"E|) 



and simulation of a larger sample of test particles, but a residual imbalance is unavoidable. 
Since the quantities we are interested in, such as frequencies and damping rates, are weak 
functions of the number of condensate atoms, a small residual drift in the condensate number 
will not affect our results significantly. 



B. Monopole modes 

This section presents the main results of the paper, where we simulate the monopole 
"breathing" mode in an isotropic trap. These calculations are not motivated by experiment 
which have yet to be performed in this geometry. Rather, we are mainly interested in 
comparing our results to previous theoretical approaches for Cu and Landau damping which 
have relied on spherical symmetry. It should be emphasized that our calculations do not 
face this restriction, though the simple geometry does allow us to more readily observe and 
quantify effects ensuing from C 2 2 and Ci 2 collisions between atoms. In fact, as reported 
elsewhere []28|-j30|1 , our methods have already been applied successfully to other experiments 
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in anisotropic traps, most notably to the study of scissors modes in which a full 3D simulation 
is necessary. 



1. Static thermal cloud approximation 

As an important test of our treatment of G\i collisions, we evaluate the damping of 
the monopole condensate mode within the so-called static thermal cloud approximation 
discussed by Williams and Griffin (WG) ||40|| . In this approximation, one considers the 
dynamics of the condensate in the presence of a static equilibrium distribution of thermal 
atoms. Due to the condensate oscillation, the condensate is no longer in local equilibrium 
with the noncondensate and as a result, Cyi collisions play a role in damping the mode. 
This effect enters through the R term in the generalized GP equation (0). It should be 
emphasized that R is provided by the theory and the relaxation it gives rise to is not 
introduced in a phenomenological way as is sometimes done [[11,32]. Linearization of the 



GP equation leads to generalized Bogoliubov equations which can be solved to determine 
collective mode frequencies and damping rates. The latter are of particular interest since 
they are directly related to the transfer of atoms between the condensate and thermal cloud 
as a result of Cyi collisions. The results obtained EDI are in fact close to those found in the 



TF approximation which gives the damping rate [H3 



hJdr6n 2 Ar)/r'(r) 

li = - 



2 <r) 



I f dr5n 2 

where Srij(r) is the density fluctuation associated with the mode j and 1/r' = gY\ 2 /k,BT. 
One sees that the damping in the TF approximation is given by a weighted average of the 
equilibrium Cyi collision rate. 

Our simulation of the static thermal cloud approximation involves the propagation of 
the condensate wavefunction according to (|3[) but with a stationary noncondensate mean 
field, 2gn (r). At the same time, the thermal atoms evolve in an effective potential defined 
by the condensate and noncondensate equilibrium densities. Although the thermal atoms 
are not allowed to undergo collisions, their dynamical evolution allows one to perform a 
Monte Carlo sampling of phase space in order to generate the Cu collision probabilities at 
each time step. These probabilites are then used to calculate the imaginary term, R(r,t), 
in the GP equation according to (|38|). These simulations can be compared directly with 



the calculations by WG KU| and therefore provide a direct test of our simulation methods, 
in particular, the calculation of Cu collision probabilities. It is important to quantify the 
errors that arise since they will also enter into our full simulations in which the effects of 
mean fields and collisions on the thermal cloud are included completely. 

The monopole mode is excited by initially scaling the equilibrium condensate wave- 
function, $(r, 0) = ct~ 3 / 2 $ ( r /ot), where the scale parameter a is 0.95. This dilation 
of the wavefunction is an alternative to imposing an initial velocity field as discussed in 



Sec. |III PL The widths of the condensate wavefunction in the x, y, and z directions are 



defined by mean-squared deviations, e.g. o x = a/ (x 2 ) — (x) 2 , where the moments are given 
by (x) = w I d r X n c( r )- Plots of these widths show a damped oscillation, and to quantify 
the frequency u> and damping rate L, we fit the data to an exponentially decaying sinusoid. 
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Since each direction gives slightly different values due to statistical fluctuations, we average 
over the three to obtain values for u and V . Our numerical results are plotted with those of 
Ref. [fHJ in Fig. f|. We find excellent agreement between the two approaches, except for the 
damping rate at T = 50 nK which is somewhat lower than the WG result. This discrepancy 
arises through errors in estimating the phase-space density in the condensate surface region 
where the fugacity approaches unity and the distribution function / is sharply peaked in 
momentum space around p = 0. The Cxi collision rate in this region is similarly enhanced, 
especially at higher temperatures. Our binning procedure is of insufficient accuracy to fully 
capture this peak, and since the surface region is the major contributor to the Cxi damping, 
this then leads to an underestimate of the rate. We illustrate this point in Fig. [| by plotting 
the result (open circle) of a simulation at T = 50 nK which uses the analytical expression 
for /o in (|39|), as opposed to the binned phase-space density. We now find much better 
agreement with the WG damping result. The generally good agreement with WG for the 
frequency and damping rate confirms that collision rates can be reliably calculated using 
our Monte Carlo sampling methods. 

Although the binning procedure introduces some minor errors into our simulations within 
the static thermal cloud approximation, we expect them to be even less important when the 
full dynamics of the thermal cloud is included. Due to mean-field interactions with the 
condensate, the thermal cloud will be strongly perturbed in the surface region and the 
distribution in phase space will tend to be "smeared out", making the binning procedure 
more reliable. C22 collisions compete against this effect by rethermalizing the particles to 
a Bose distribution; however, this can only make a significant difference if the collisional 
timescale is short compared to that of the oscillation. For the present calculations, we have 
^0^22 >> 1 and the gas is in the collisionless regime. We would therefore expect the thermal 
cloud dynamics to be very important in determining the damping due to Cxi collisions. 

To illustrate this we have performed full simulations including mean-field interactions and 
collisions at T = 20 nK and 30 nK. The results obtained with only C22 collisions included 
are shown by open squares, while the results including Cxi collisions as well are shown by 
the full squares. One sees that the overall damping rate increases by only 5-10% when Cxi 
collisions are added in. In fact, collisions of either kind contribute little to the damping which 
is dominated by Landau damping (as discussed in the following subsection). Furthermore, 
we find a small downward shift in the frequency compared to the zero temperature value, 
in contrast to the significant increase seen within the static approximation. This increase is 
due to the fact that the condensate is oscillating in the presence of the static mean field of 
the equilibrium thermal cloud which effectively enhances the oscillator frequency of the trap. 
This effect is eliminated when the thermal cloud is allowed to respond to the dynamic mean 
field of the condensate. As we shall further demonstrate in the following subsection, dynamic 
mean-field effects typically dominate the finite temperature behavior, with collisions playing 
a secondary but important supporting role in equilibrating the system. 

As regards the status of the static thermal cloud approximation, we conclude that it pro- 
vides a useful method for qualitatively determining the effects of Cxi collisions on collective 
modes. However, its quantitative predictions for mode frequencies and damping rates are 
unreliable. 
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2. Landau damping 



As our final example, we have performed simulations for the system studied by Guilleu- 
mas and Pitaevskii P4j] , namely 87 Rb atoms confined in an isotropic trap of frequency 
ujq = 2tt x 187 Hz. To begin, we consider a total of TV = 5 x 10 4 atoms and excite the 
monopole mode by an initial scaling of the condensate radius by a factor of a = 0.9, with 
the thermal cloud initially in its equilibrium state. The condensate width oscillations are 
then followed over a time scale of ou^t = 30. Fig. |5] shows damping rates and frequencies 
as a function of temperature found by fitting an exponentially decaying sinusoid to the 
time-dependent width. At each temperature three simulations are performed. The first in- 
volves free propagation of thermal test particles without collisions, corresponding to solving 
the collisionless Boltzmann equation. The second includes C22 collisions between thermal 
atoms, while the third includes both C22 and C12 collisions. At low temperatures, all three 
simulations give similar results, reflecting the fact that the number of thermal atoms is small 
and collisions play a minor role. With increasing temperature, the differences between the 
simulations increase. Qualitatively, the behaviour is similar to what was found previously 
for the scissors mode p9| ; collisions have the effect of shifting the frequency downward as 
compared to the collisionless result, and significantly enhance the damping rate. The effect 
of C22 collisions is particularly strong at high temperatures, which at first sight may seem 
surprising since C22 collisions do not couple to the condensate directly. 

To gain more insight into the collisional dependence of the damping, we focus on the 
time-dependent evolution for a particular temperature, T = 200 nK (T/T 6 ° = 0.644, where 

1 /3 

T c ° = 0.94/j£ t ;oA tot [||] is the transition temperature of the corresponding ideal gas in the 
thermodynamic limit). Fig. ^ plots a x vs. t for the collisionless and full (C12+C22) simula- 
tions. The initial damping rate in both calculations is seen to be similar, however at later 
times the collisionless oscillation departs from a simple exponential decay, and the oscillation 
amplitude tends to saturate. This behavior is not seen to the same degree when collisions 
are included. To quantify this behavior, we define a local damping rate by fitting a damped 
sinusoid to the data within a window of width A(u t) = 9 centered on the time t. Fig. [7| 
plots this local damping rate as a function of t. We see large variations in the damping rate 
with time, with the largest rate occuring initially. The deviations from the initial value are 
largest in the collisionless case, where the damping rate dips nearly to zero. Similar behavior 
is observed over the whole range of temperatures, and accounts for the lower damping rates 
obtained by fitting the entire data set. 

To explain this behavior, we note that damping of the condensate oscillation is associated 
with the transfer of energy from the condensate to the thermal cloud. If this energy exchange 
is mediated by mean-field interactions, it is referred to as Landau damping. From the point 
of view of the thermal cloud, the dynamic condensate mean field 2gn c (r, t) acts as an external 
perturbation which can lead to the excitation of thermal atoms. Of course, the rate at which 
these excitations occur depends on the phase-space distribution of the thermal particles. In 
our simulations, the thermal cloud is initially in an equilibrium state and the damping 
rate is observed to be independent of collisions. This damping is essentially pure Landau 
damping and its magnitude is determined by the rate at which the oscillating condensate 
can do work on the equilibrium thermal distribution. In this respect, our initial damping 
rate is analogous to conventional perturbation theory estimates (as discussed below). As 
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time progresses in our simulations, the thermal cloud begins to deviate from an equilibrium 
distribution and the magnitude of Landau damping is correspondingly affected. Evidently, 
the perturbation of the thermal distribution is such as to reduce the rate of energy transfer 
to the thermal atoms, whereupon the damping rate decreases with time as seen in Fig. [7]. 
The deviation is in fact a nonlinear effect as it was found to depend on the amplitude of the 
condensate oscillation. An analogous effect appears in the context of plasma oscillations, 
where Landau damping is due to the energy transfer from the collective plasma wave to 



single-electron excitations [45.46 



In the absence of collisions, the distribution of thermal atoms continues to evolve in a 
complicated way and the effective damping rate exhibits an oscillatory time dependence. 
However, as soon as C22 collisions are switched on, the damping rate deviates less strongly 
from its initial value. The effect of these collisions is to drive the thermal cloud towards a 
state of local equilibrium and the damping rate tends to maintain its original value. The 
inclusion of C12 collisions has a similar effect and we find a damping rate which is almost 
time-independent when both collision processes are retained. However, Cxi collisions do more 
than simply equilibrate the thermal cloud since they also lead to the source term R(r, t) in 
the GP equation. As we have already discussed, this term gives rise to its own contribution 
to damping which is quite separate from Landau damping. It should be emphasized that 
it is impossible to separate the total damping rate into individual components. Mean-field 
and collisional effects are interrelated, and all must be included to completely account for 
the actual damping rates. 



We next turn to a comparison of our results with those of Guilleumas and Pitaevskii [44 



Since we have used quite different methods to calculate damping rates, it is useful to first 



discuss the perturbation theory calculation of Landau damping used by these authors [47 



Within this approach, Landau damping is associated with the decay of a mode of oscillation 
(with energy hou osc ) as a result of the excitation of a thermal quaisparticle from an initial 
state of energy Ei to a final state of energy E k . The damping rate is then given by Fermi's 
golden rule ||,0,||,|T|] 



-y 

ik 



Aik\ [f(Ei) - f(E k )}5(E k -Ei- fuj osc ), (41) 



where the sum is over all excitations that satisfy energy conservation, while the matrix 
element A ik depends upon the form of the excitations. The thermal states are occupied ac- 
cording to the equilibrium Bose distribution f(E). This damping rate is therefore analogous 
to the initial damping rate we obtain in simulations which start with an equilibrium thermal 
distribution. 

Guilleumas and Pitaevskii ||4| evaluate ((41] ) as a function of temperature using Bogoli- 
ubov excitations of the condensate as the thermal quasiparticles. These are determined from 
the Bogoliubov equations for a fixed number of condensate atoms, N c ; the corresponding 
number of thermal atoms is then a function of temperature and is given by summing over 
the thermal occupation f(Ei) of the quasiparticle states. To actually evaluate the Landau 
damping rate at the frequency of the monopole mode of interest, the delta functions in ( f4T|) 
are replaced by Lorentzians of width A. They show that the results obtained are essentially 
independent of this parameter. 

To compare with these results we performed collisionless simulations and extracted the 
initial damping rate as discussed earlier. The comparison is made in Fig. |] where results are 
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presented as a function of temperature for N c = 5 x 10 4 and 1.5 x 10 5 . Given the completely 
different methods of calculation, the agreement is remarkable. The agreement persists even 
down to low temperatures where one might expect differences to appear as a result of our use 
of semiclassical HF excitations as opposed to the Bogoliubov excitation spectrum. The fact 
that the agreement is as good as it is is perhaps understandable in view of the observation 
in Ref. p8| , |49| that the density of states in the HF and Bogoliubov approximations are very 
similar. Although the semiclassical HF approximation was not discussed, we would of course 
expect it to be close to the quantal HF result. Since the density of thermal excitations is 
an important ingredient in the calculation of Landau damping, we can begin to see why our 
semiclassical calculations give very similar results. 



V. CONCLUSIONS 

In this paper we have provided a detailed description of the numerical scheme we have 
used to simulate trapped Bose-condensed gases at finite temperatures, based on the ZGN for- 
malism which treats the thermal excitations semiclassically within a Hartree-Fock-Popov ap- 
proximation. The procedure involves solving simultaneously a generalized Gross-Pitaevskii 
equation for the condensate and a Boltzmann kinetic equation for the thermal cloud. The 
two equations are coupled by mean fields and collisions, both of which influence the dy- 
namics of the two components in significant ways. Our scheme has been carefully tested to 
ensure that it provides an accurate description of the system dynamics. In particular, we 
have shown that iV-body simulations, together with the Monte Carlo sampling of collisions, 
is an effective and reliable method for determining the thermal cloud dynamics. 

Our scheme can be used to model the dynamics of the gas over a wide range of temper- 
atures and physical conditions. As an example, we have studied the monopole "breathing" 
mode in a spherical trap. Two sets of calculations were performed. The first provided a check 
of the treatment of collisions within the static thermal cloud approximation of Williams and 



Griffin [iO]. Results within this approximation were reproduced, but full dynamical sim- 
ulations indicated that the approximation is primarily useful as a qualitative indicator of 
the effect of C\2 collisions. Unfortunately, its quantitative predictions for mode frequencies 
and damping rates cannot be trusted. Our second set of simulations focussed on Landau 
damping. This is typically the dominant damping process for condensate modes at finite 
temperatures. However, the damping observed in a simulation, and by extension in real 
experimental situations, is determined by a delicate interplay of the mean-field excitation 
of the thermal cloud and collisions. The thermalizing effect of the latter strongly influences 
the rate at which mean-field excitations take place. 

We also compared our results for Landau damping to those of Guilleumas and Pitaevskii 



44j, and very good agreement was found. This confirms that the semiclassical HF descrip- 
tion of the thermal cloud reproduces the Landau damping as calculated using Bogoliubov 
excitations. This is not too surprising since the density of excitations in the two approx- 
imations is very similar. However, as we have already explained, the Landau damping as 
determined by assuming the thermal cloud to be in an equilibrium state is not necessarily 
the damping that will be observed in an experiment. A consistent treatment of the dynamics 
of the condensate and thermal cloud is needed in order to make detailed comparisons with 
experiment. 
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Applications of our technique to other collective modes have also been made and are 
discussed elsewhere [28,P9| . Our results for the temperature-dependent damping and fre- 



quency shifts are in good agreement with experiment for both quadrupole J28J and scissors 
p9| modes. This in itself confirms the accuracy of our theoretical formulation of the system 
dynamics and the numerical methods used. It is hoped, however, that the present paper 
provides more insight into the content of the theory and the reasons for its success. Interest- 
ing future systems for study could include topological defects (e.g. vortices and skyrmions), 
optical lattices, and dynamical instabilities of surface modes in the presence of a rotating 
thermal cloud pD[ . 

We conclude with a few comments about where we may go next. It would be useful 



to incorporate Bogoliubov excitations in the kinetic theory in place of HF excitations [pi 



Although this does not seem to be important for Landau damping, the HF approximation 
may not be good for certain situations where the thermal occupation of the lowest modes 
are significant (e.g., for large atom numbers or at low temperatures). Another possibility 
is to implement a hybrid scheme, where highly-occupied, low-lying modes are treated using 
classical field methods |5^ , |53| , while the rest are treated semiclassically using the present 



technique. Finally, it would be of interest to investigate the importance of going beyond 
the Popov approximation by including the anomalous average, m. In doing so one must be 
careful to ensure that the new model is gapless 0. This may involve renormalization of the 
coupling constant g |Q or replacing the contact potential by a generalized pseudopotential 
E9. 
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APPENDIX A: MONTE CARLO CALCULATION OF COLLISION RATES 

Our purpose here is to show how a Monte Carlo evaluation of the collision rate in (|28D 
leads to a definition of collision probabilities to be used in the simulations. We first note 
that d 3 rd 3 p/h 3 C%2 t represents the number of atoms leaving the phase space volume element 
d 3 rd 3 p/h 3 per unit time as a result of collisions. Integrating this over momenta gives the 
number of atoms in d 3 r suffering a collision per unit time. Thus the mean collision rate per 
atom and per unit volume is 

rS* = / ^cif , (Ai) 

which is the quantity displayed in (^7]). We now write the required local collision rate as 

= J d 6 pw(p)g(p) (A2) 
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where p is a point in 6-dimensional momentum space and the factor w(p) = / '(pi) / (P2) / h 6 is 
considered as a weight function. We denote the maximum value of w(p) by w max and define 
the domain on which the integrand is nonzero by [— p max /2, £> m ax/2] for each momentum 
component. Choosing a point Pi at random in the hypervolume (p max ) 6 , and a random 
number Ri uniformly distributed on [0, to ma x], the point Pi is accepted if Ri < w(pi) and the 
quantity g(pi) is accumulated. The value of the integral is then given approximately as 

T°If ~ (p max ) 6 w max ^J^g( Pl ) , (A3) 

i 

where A" is the number of random pi points chosen and the prime on the summation includes 
only those points for which Ri < w(pi) . For g = 1, the integral is simply n(r) 2 . Thus, 

n{rf = (p max ) 6 «W^ (A4) 
where N s is the total number of points accepted, and 

r " * ^jfU'gipi) • (as) 

i 

The sample of A" s points accepted consists of N s pi-values and A^ p 2 -values, each of 
which is distributed according to /(p). This set of 2N S p- values can be identified with Af ce ii 
test particles in a cell of volume A 3 r. If this set is to be representative of the local density, 
we must have 

iV cell 2N S 

" (r) = = • (A6) 

With this identification, 

A 3 rT^ = 2nJ29(PlP2)- (A7) 

8=1 

In other words, the collision rate can be estimated by sampling the test particles in the cell 
A 3 r in pairs. Inserting the explicit form of g for the 22 collision rate in (p8|), we have 

A 3 rr- t = 2^n(r) ( x|v i -v,| f ^(1 + / s )(l + / 4 ) , (A8) 

where the sum is now taken over pairs of test particles. This expression allows us to define 
the probability P? 2 that a pair of atoms (ij) in the cell suffers a collision in a time interval 
At: 

P% = nOMv, - v,| J ^(1 + / s )(l + / 4 )At . (A9) 

Selecting atoms in pairs from each cell and assigning them a collision probability P? 2 allows 
us to simulate the effect of collisions in a way which is consistent with the Boltzmann collision 
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integral. Note that the factor of 2 in ( |A8|) accounts for the fact that 2 atoms are affected for 
each pair collision. This factor is therefore not included in the definition of the pair collision 
probability. 

We treat C*i 2 collisions somewhat differently. First, we note that the total rate of change 
of the number of thermal atoms per unit volume due to these collisions is 

-TT^u = ^— / dp 2 dp 3 dp 4 5(mv c + p 2 - p 3 - p 4 )5(e c + e 2 - e 3 - e 4 ) 

hr 7im z h 6 J 

x[/ 2 (l+/ 3 )(l+/4)-(l + / 2 )/3/4] 

= - i% . (Aio) 



According to this definition 

7cm 2 h 3 



(jfi 

^12 = / dp 2 dp 3 dp 4 5(mv c + p 2 - p 3 - P4 )5(e c + e 2 - e 3 - e 4 )/ 2 (l + / 3 )(1 + U) 



= J d ^hn c avT j §(l + h)(l + h) (All) 

is the rate of decrease of the number of condensate atoms per unit volume as a result of a 
collision with a thermal atom (hence the designation 'out'). This rate can be estimated by 
writing 

T°f = j d 3 p 2 w(p)g(p) , (A12) 

where w(p) = f{p)/h 3 and g(p) is the remaining part of the integrand. A Monte Carlo 
sampling of the integral leads to the estimate 



A 3 rr?f^^), (A13) 



i=i 



where N s represents the number of atoms in the cell of volume A 3 r. The probability of an 
atom in the cell suffering this kind of collision in the time interval At is therefore 

Pr = g(pi)At = r^xyVc ~ v^| 2 - Agnjm J ^(1 + / 3 )(1 + / 4 )At , (A14) 



which is the origin of the expression given in (35) 
The 'in' collision rate is given by 



r i n 2 = — vr^ I dp 2 dp 3 dp 4 5(mv c + p 2 - p 3 - p 4 )5(e c + e 2 - e 3 - e 4 )(l + / 2 )/ 3 / 4 
itm z h 6 J 

-firh / ~h3~f*^ m ~ s ^ Pc ~ p ^ ' ( pc _ P2 ^ ~~ m 9 n c](± + is) , (A15) 

where we have interchanged the particle labels 2 and 3 to obtain the second line in this 
equation. This rate corresponds to two thermal atoms scattering into a condensate atom 
and an outgoing thermal atom, and is thus the rate that atoms feed into the condensate as 
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a result of collisions. Although the collision of atoms 2 and 4 can be treated by the methods 
used to analyze the C22 collision rate, it is preferable to define a single atom collision rate 
by writing this integral in the form of (|A12[) and performing a Monte Carlo sampling with 
respect to the p 2 variable. This procedure leads to the collision probability per atom 

Pr = At f + f 3)h T ^ 6 [(p c - p 4 ) • (p c - p|) - mgn c } , (A16) 

which is simplified and discussed further in the body of the paper. 
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FIGURES 
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FIG. 1. Equilibrium noncondensate density against position, along a line through the center of 
an isotropic trap with frequency ujq = 2ir x 187 Hz. The system consists of a total of iV-tot = 5 x 10 
87 Rb atoms at a temperature of T = 250 nK. The critical temperature for an equivalent ideal 
gas would be T® = 310.6 nK. Nj> = 4.0 x 10 5 test particles are sampled according to the actual 
equilibrium density (solid line) . The fluctuating dashed line is a result of binning particles using a 
cloud-in-cell method, while the smooth dashed line shows the effect of convolving the cloud-in-cell 
density with a Gaussian. 
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FIG. 2. The T Q 22 collision rate as a function of position, for an equilibrium distribution and the 
same parameters as Fig. |]. The solid line shows the result of a direct evaluation of (p9|), while the 
points plot the results of a Monte Carlo evaluation (|3C|), 
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FIG. 3. Same parameters as in Fig. [I], for T 12 collisions between the condensate and thermal 
cloud in equilibrium. The solid line plots a direct evaluation of (|34|), while the circles shows a 
Monte Carlo calculation for the "out" rate (solid) and "in" rate (open). 
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FIG. 4. Temperature dependent (a) frequency shifts and (b) damping rates of the condensate 
monopole mode in a spherical trap (ujq = 2tt x 10 Hz), in the presence of a static thermal cloud. 
The total number of atoms is TV-tot = 2 x 10 6 . The critical temperature for the corresponding ideal 
gas is Tq = 56.8 nK. Our results are plotted as solid circles, while the solid line is the prediction 
of Williams and Griffin |^(|. The open circle at T = 50 nK is the result of a calculation using the 
analytical form for the phase-space density (p9|). The squares plot results of simulations including 
thermal cloud dynamics, with C22 collisions only (open) and both C22 and C12 collisions (closed). 
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FIG. 5. (a) Frequency and (b) damping rate of a monopole mode (o>o = 2ir x 187 Hz, 
-^tot = 5 x 10 4 ), including thermal cloud dynamics. Results are shown for simulations with no 
collisions, C22 = C12 = (closed circles), C22 collisions only (open circles), and both C12 and C22 
collisions (inverted triangles). 
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FIG. 6. Time-dependent width of the condensate, a x , after excitation of the monopole mode at 
T = 200 nK. The dashed line shows the collisionless evolution, while the result of a full simulation 
(C12 and C22) is indicated by the solid line. 
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FIG. 7. Damping rates for the same parameters as in Fig. |(| where fits are taken with a series 
of windows in the range [u>ot — 4.5, u>ot + 4.5]. We plot data for simulations which include no 
collisions (open circles), Cvi only (inverted triangles), C22 only (solid triangles), and both C12 and 
C22 collisions (solid circles). 



3 





at 


0.06 - 






V' 


0.04 - 


*/ / 






0.02 - 








0.00 - 


1 1 1 1 



0.0 0.5 1.0 1.5 2.0 2.5 

kT//LL 

FIG. 8. Initial damping rate calculated over the first time interval of size A(uot) = 9 
(points) compared to the results of Guilleumas and Pitaevskii |44j (lines). Results are plotted 
for N c = 5 X 10 4 (solid points and line) and -/V c = 1.5 x 10 5 (open points, dashed line) condensate 
atoms. Following Ref. p4j quantities are plotted in terms of dimensionless units, with the ratio 
T/ujm (damping rate over mode frequency) plotted against hsT/fi. For T/ujm we calculate the 
mean over the three directions, while the standard deviation yields a rough estimate of the error. 
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